Creative Assignment 4

library(osmdata)
library(opentripplanner)
library(tidyverse)
library(sf)
library(ggthemes)
library(ggspatial)

Load locations for emergency shelters

For this assignment, I am interested in looking at the walksheds and drivesheds to Emergency Shelters in Cambridge. It is important for Cambridge residents to be able to quickly access these shelters during emergency situations.

Emergency Shelters data from: https://www.cambridgema.gov/GIS/gisdatadictionary/Public_Safety/PUBLICSAFETY_EmergencyShelters

Reading the shapefile that was saved to my computer:

emergency_shelters <- st_read(
  "data_downloads/PUBLICSAFETY_EmergencyShelters.shp/PUBLICSAFETY_EmergencyShelters.shp")

## Reading layer `PUBLICSAFETY_EmergencyShelters' from data source `/Users/jocelyntsai/Desktop/GitHub/ctsai1-vis/data_downloads/PUBLICSAFETY_EmergencyShelters.shp/PUBLICSAFETY_EmergencyShelters.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15 features and 12 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 749626.1 ymin: 2954843 xmax: 767879 ymax: 2968207
## projected CRS:  NAD83 / Massachusetts Mainland (ftUS)

The data that I downloaded was originally in NAD 1983 StatePlane Massachusetts Mainland coordinate system, I have to transform it to latitude/longitude first using WGS84 (EPSG:4326) bounds or else I will run into en error when trying to make isochrones later.

emergency_shelters <- emergency_shelters %>%
  st_transform(crs=4326)

Get street network data

The opq function creates a call to the open street map server to download street network data. bbox (bounding box), is like the search bar on google maps (typing in general area name is fine, it doesn’t have to be exact). osmdata_xml will create a file with the information about the street network, and it will save it to the default directory created.

opq(bbox = 'Cambridge MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/cambridge_streets.osm')

I am transforming my spatial data of street lines onto the NAD83(HARN) / Massachusetts Mainland projected coordinate system from https://spatialreference.org/ref/epsg/2805/ . The add_osm_feature when ‘highway’ is specified takes all the roads in Open Street Map. I then used osmdata_sf to also download data in the sf simple features format so I can plot it.

MA_mainland <- "+proj=lcc +lat_1=42.68333333333333 +lat_2=41.71666666666667 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

cambridge_street_features <- opq(bbox = 'Cambridge MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

cambridge_streets <- cambridge_street_features$osm_lines %>%
  st_transform(crs = MA_mainland)

I then plotted the streets in Cambridge to see what’s included:

ggplot(cambridge_streets) +
  geom_sf() 

Set up Open Trip Planner

Download a Java utility called otp.jar and save it to the OPT folder (only need to run this code chunk once):

path_otp <- otp_dl_jar("OTP")

## The OTP will be saved to OTP/otp.jar

Build a graph representing street and transit networks:

path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024) 

Setting up Open Trip Planner (OTP application opens in the web browser):

otp_setup(otp = path_otp, dir = path_data, memory =1024)

Connect to Open Trip Planner:

otpcon <- otp_connect()

## Router http://localhost:8080/otp/routers/default exists

Create isochrones

I created isochrones for areas within an 8 minute walk and an 8 minute drive of emergency shelters in Cambridge using the otp_isochrone function. I then used the rbind function to combine these sets of polygons (drive, walk) into one dataframe. otp_stop() function stops the Open Trip Planner application from running.

iso_8min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters, 
                mode = "WALK", cutoffSec = 480) %>%
  st_transform(crs = MA_mainland) %>%
  mutate(mode = "walk")

iso_8min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters, 
                mode = "CAR", cutoffSec = 480) %>%
  st_transform(crs = MA_mainland) %>%
  mutate(mode = "drive")

## Warning in otp_isochrone(otpcon = otpcon, fromPlace = emergency_shelters, :
## Failed to get isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7ff1"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fee"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fe9"}]}Failed to get
## isochrone with error: {"type":"FeatureCollection","features":
## [{"type":"Feature","geometry":null,"properties":{"time":
## 480},"id":"fid--6db43113_17501b2ddf3_-7fe8"}]}

iso_all_modes <- rbind(iso_8min_walk, iso_8min_drive)

otp_stop()

Map isochrones

Draw the walk/drive isochrones using the street network as a basemap

left_side <- st_bbox(iso_all_modes)$xmin
right_side <- st_bbox(iso_all_modes)$xmax
bottom_side <- st_bbox(iso_all_modes)$ymin
top_side <- st_bbox(iso_all_modes)$ymax

ggplot(iso_all_modes) +
  geom_sf(data = cambridge_streets, color = "gray") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "plasma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map()

Draw the walk/drive isochrones using OpenStreetMap as a basemap

Typing rosm::osm.types() into the console shows all other basemap options.

1. default

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "magma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

## Loading required namespace: raster

2. hillshade

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "hillshade", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "plasma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

3. cartodark

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "cartodark", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking")) +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

4. cartolight

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "plasma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

5. thunderforestoutdoors

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "thunderforestoutdoors", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "plasma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

6. osmgrayscale

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "osmgrayscale", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking"),
                       option = "plasma") +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

7. thunderforestlandscape

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type= "thunderforestlandscape", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.3) +
  geom_sf(data = emergency_shelters) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 8 minutes",
                       labels = c("By Driving", "By Walking")) +
  annotation_scale(location = "br")+
  annotation_north_arrow(location = "tr", style= north_arrow_minimal())+
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Isochrone areas: scatter plot

To calculate the area of each isochrone, I used the st_area() function. The pivot_wider() function then creates a seperate column for each mode of transportation (drive or walk), and each row represents a location with both the drive and walk isochrones instead of having each row in the dataframe just containing either one of the isochrones. For labels, I used breaks/1000000 to convert meters square to kilometer square.

iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(drive))) +
  geom_point() +
  scale_x_continuous(name = 
            "Area within an 8 minute walking distance\nof an emergency shelter (square km)",
            breaks = breaks <- seq(10000, 500000, by = 50000),
            labels = breaks / 1000000) +
  scale_y_continuous(name = 
            "Area within an 8 minute driving distance\nof an emergency shelter (square km)",
            breaks = breaks <- seq(0, 3000000, by = 500000),
            labels = breaks / 1000000)+
  theme_economist()

## Warning: Removed 4 rows containing missing values (geom_point).

Isochrone areas: scatter plot with polar coordinate system

ggplot(iso_areas, 
       aes(x = as.numeric(walk),
           y = as.numeric(drive))) +
  geom_point(alpha = 0.8, size = 1) +
  stat_smooth(color = "lightblue", linetype = 1, size = 0.8) +
  scale_x_continuous(name = "Area within an 8 minute walking distance\nof an emergency shelter (square km)",
                     breaks = breaks <- seq(10000, 500000, by = 50000),
                    labels = breaks / 1000000) +
  scale_y_continuous(name = "Area within an 8 minute driving distance\nof an emergency shelter (square km)",
                     breaks = breaks <- seq(0, 3000000, by = 500000),
            labels = breaks / 1000000)+
  coord_polar(theta = "x") +
  theme_minimal() 

## Warning: Removed 4 rows containing non-finite values (stat_smooth).

## Warning: Removed 4 rows containing missing values (geom_point).